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Abstract. We studied the pure and dilute Baxter- Wu (BW) models using the Wang- 
Landau (WL) sampling method to calculate the Density-Of-States (DOS). We first 
used the exact result for the DOS of the Ising model to test our code. Then we 
calculated the DOS of the dilute Ising model to obtain a phase diagram, in good 
agreement with previous studies. We calculated the energy distribution, together with 
its first, second and fourth moments, to give the specific heat and the energy fourth 
order cumulant, better known as the Binder parameter, for the pure BW model. For 
small samples, the energy distribution displayed a doubly peaked shape, and finite size 
scaling analysis showed the expected reciprocal scaling of the positions of the peaks 
with L. The energy distribution yielded the expected BW a = 2/3 critical exponent for 
the specific heat. The Binder parameter minimum appeared to scale with lattice size 
L with an exponent 9b equal to the specific heat exponent. Its location (temperature) 
showed a large correction-to-scaling term 6\ = 0.248 ±0.025. For the dilute BW model 
we found a clear crossover to a single peak in the energy distribution even for small 
sizes and the expected a = was recovered. 
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1. Introduction 

The two-dimensional Ising model has received such widespread attention as the 
paradigm system for phase transitions, that one sometimes says that a certain system 
is the "Ising model" of a class of problems. While it is clearly special, it is not the only 
two dimensional model of phase transitions with an exact expression for its free energy 
Its critical behavior is, in fact, rather atypical relative to many other two-dimensional 
systems and even to the three-dimensional Ising model, especially in the specific heat, 
where its critical exponent, a, is zero. The nature of the corrections-to-scaling in the 
spin 1/2 Ising model is also very different to that of many other interesting systems. 

Another spin system, now known as the Baxter Wu (BW) model, was solved by 
R.J. Baxter and F.Y. Wu [HI21- Spins (jj = ±1, are situated on the triangular lattice 
and interact via a three spin interaction, 

H = -J^CTiCTjO-k, (1) 

where i,j and k are the vertices of a triangle as shown in Fig. ^ J > is the 
ferromagnetic coupling between nearest neighbor spins. The BW model exhibits a 




Figure 1. The energy of a given configuration is the sum of all interacting triangles 
formed by nearest neighbor spins 

second order phase transition with its critical temperature (T c ) given by 2J/kT c = 
ln(l + V2) = 2.26918..., (the same numerical value as for the Ising model on the square 
lattice). The specific heat critical exponent is equal to the correlation length exponent, 
a — v — 2/3. Series-expansion results gave the conjectured magnetization exponent 
of /3 = 1/12 and a susceptibility exponent of 7 ~ 1.17 [4J. The latter confirmed the 
prediction of 7 = 7/6 from the well known scaling relation a + 2/3 + 7 = 2 (SlIEl- 
Real Space Renormalization Group methods have also been used [3 |Hl E] to study 
the pure model, and the critical eigenvalues obtained gave critical exponents consistent 
with series-expansion and exact results. An exact form for BW corrections-to-scaling 
was found by Joyce JU] who conjectured that the spontaneous magnetization varied 
as M = t@ (fo(t) +t 2 ^fi{t) + • • •) with analytic functions fo, fi of the distance t = 
(T — T c )/T c . Adler and Stauffer confirmed this with series and Metropolis Monte Carlo 
estimates jTTj. 
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Dilute Ising models are also somewhat famous, but for rather different reasons, 
as they have been the source of a great deal of controversy. Presumably because of 
the anomalous specific heat structure in the pure case, numerical work in the dilute 
regime, especially near the pure limit is painful, and although a majority of authors 
(see e.g. Roder et al, [E]) have claimed that the controversy is resolved in favor of SSL 
theory [T31 EH Hi] , more study is useful. 

The annealed dilute BW model was studied by Kinzel, Domany and Aharonv[To] 
who showed from this exploration that its dominant critical behavior is in the 
universality class of the four state Potts model, although the Potts model has logarithmic 
correction terms for this case. Domany and Riedel |17j . argued the same for the pure 
BW model by means of symmetries of the Landau- Ginzburg- Wilson Hamiltonian. (Note 
that there are first order fixed points in the neighborhood of these models). The 
quenched dilute BW model was studied by Landau and Novotny [THj, who found a 
substantial change in the critical behavior of the specific heat ^H] for an impurity 
concentration ofl — x = 0.1. They also conjectured that the zero temperature threshold 
concentration above which no long-range order could be seen was about x c ~ 0.71. 
(See also results from a cluster-algorithm study of this system ;20j). More recent 
calculations [2H 123 showed that the value of x c is even higher (x c ~ 0.755), and is 
bounded by x l ° w = 0.710 ± 0.001 and x^ igh = 0.784 ± 0.004. This is substantially 
above the value for Ising models where x c is simply the percolation threshold of the 
corresponding lattice, which is rarely above 0.5 . 

Recently Wang and Landau E3j proposed a very efficient algorithm for 
calculating the density-of-states (DOS), (i.e. the degeneracy of any level in energy 
space), g(E), for Ising models and some related systems. To explore the issues of both 
pure and dilute BW models further, and to see how its different particulars of large a 
and corrections to scaling emerge from the calculation of the DOS, we have chosen to 
apply the WL algorithm to both the pure and quenched dilute BW models, and to study 
the behavior of the energy distribution and related moments [23 I2H| using the simulated 
DOS. The DOS of the pure and dilute Ising models was studied for comparison purposes. 

In the next section we discuss the WL algorithm. In section E] we present a 
comparison of an exact calculation of the DOS for the Ising model j22j with simulations 
using WL and give some results for the dilute Ising case. In section HI we give in detail 
our results for the pure BW model, and in section El the results for the dilute BW model 
are presented. Finally we discuss the implications of our results in section El 

2. The simulation method 

Conventional Monte-Carlo (MC) methods [2B1 EH ED] generate the canonical energy 
distribution at a given temperature T . It is usually narrowly peaked around 
this temperature. The need to perform multiple simulations in order to obtain 
thermodynamics in a large range of temperatures requires a large computational effort. 
Other methods based on histogram accumulation [213 EI] approximate the distribution 
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by the energy histogram at T . This distribution can then be reweighted to give statistics 
at another temperature. The reweighted distributions, however, are also restricted to a 
very narrow range of temperatures and suffer from large statistical errors in their tails 
for temperatures far from T . The broad histogram method |32j calculates the DOS 
through the consideration of the average number of visits to any two adjacent energy 
levels. Lee offered the entropic sampling method using the observation that if the 
transition probability between any two energy levels is proportional to the ratio between 
the DOS of these levels, then a crude estimate to the DOS can be given when sampling 
at infinite temperature. 

Wang and Landau improved Lee's method by introducing a modification factor 
which together with generating a "flat" histogram (we have used the condition \H(E) — 
(H)\/(H) < 0.05 for any E), carefully controls the updating of the DOS. By dividing 
the energy space into different segments, and performing an independent random walk 
in each segment, one can generate very accurately, in a reasonable amount of CPU time, 
the DOS of the whole energy space, thus obtaining the canonical distribution at any 
desired temperature. 

3. Ising model results 

3.1. The pure Ising model 

We began by validating the accuracy of our implementation of the WL algorithm against 
exact results for the Ising model on the square lattice with no impurities. A detailed 
comparison was made for the case of L = 32. The partition function for the Ising model 
on a lattice of length L can be written as a low temperature expansion 



where iV = L x L is the number of spins, K = J/kT is the reduced inverse temperature, 
and x = e~ 2K is the low temperature variable. Each energy level can be labeled (relative 
to the ground state energy -2JN) by E e = Ut (£ = 0, 2, 3, ...N — 2, N), so that g t is 
its corresponding DOS. Beale [27] used an extension of Onsager's solution [34 , to give 
the exact expression for the partition function on a finite lattice |25|, and extracted 
the DOS coefficients from the expansion (|2*|) . When we plotted our results on top of 
Beale's expression for the case of L = 32 we saw no deviations between the exact and 
the simulated data within the resolution of the figure. The relative error between the 
exact and simulated data was also plotted and was found to be three orders of magnitude 
smaller than the calculated DOS and two orders of magnitude larger from the systematic 
error due to the choice of the final modification factor fa nai \ = 0.001. This showed that 
the choice of this quite large fa na \ was sufficient, so that only a relatively small number 
of iterations was required for all the simulations performed throughout this work. 

Further results from the pure Ising simulations will be introduced for comparison 
purposes in section HI 




(2) 
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3.2. The dilute Ising model 

We continued the validation process by studying the dilute Ising model, at a lattice size 
of L = 22. The Hamiltonian for the dilute Ising model may be written as 



where the random disorder variables take the values and 1, such that their 
configurational average is equal to a dilution of < x < 1. We considered the position 
of the specific heat maxima, Tc max , for different nominal concentrations centered around 
the values x = 0.8, 0, 9 and x = 0.95, as indicated in Fig. HJ It is clearly seen that for 
large concentrations (x > 0.9) the circles tend to a continuously critical line, slightly 
shifted from the solid line. The shift is a finite size effect due to the use of a small sample. 
For smaller concentrations there is a large dispersion of the circles and the data are less 
reliable. As shown in Fig. El the specific heat maximum becomes broader with decreasing 
concentration and is hard to locate precisely. The reason for this is that when lowering 
the concentration isolated clusters which rarely interact with each are formed, and hence 
energy fluctuations become smaller. For a concentration of x = 0.75 these fluctuations 
are also nearly constant and therefore no pronounced peak can be identified. It should 
be noted that presumably, when much larger samples would be used, a pronounced peak 
should be clearly seen for concentration even lower than x = 0.75 (see, for example |36j). 
In the absence of analytic results, the location of our points close to earlier estimates 
validates both our dilute code and our analysis methods. 

4. The pure Baxter- Wu model 

We calculated the DOS for the BW model using lattice sizes L ranging from 6 to 120, 
with periodic boundary conditions being imposed. For each lattice size the data was 
collected separately for each energy segment and then was combined to give the density 
of states for the entire energy landscape. We averaged over nine different runs for 
L = 30, and saw that the fluctuations were three orders of magnitude smaller then the 
measured quantity (In g), so that we neglected these fluctuations and for each lattice size 
we executed a single run per segment only. By symmetry, for any state with negative 
energy, there exist a state with positive energy, so that it was sufficient to carry out 
the random walk only for non positive energies. (A similar argument holds for the Ising 
model). Plots of the internal energy, specific heat, free energy and entropy are given in 



Early simulations jTHj showed the formation and motion of domains around the 

ferromagnetic and ferrimagnetic ground states, due to the special connectivity of the 
BW model, causing low frequency large energy fluctuations. These fluctuations made 
the impression that the system was in a metastable state, thus indicating a first order 
transition. In Fig we examined the energy distribution at Tc max and found a doubly 
peaked curve (see ref. [HH]). The system appears to fluctuate between these two 




(3) 



Fig.H 
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Figure 2. Critical line T c (x) in the T-x plane of the dilute Ising model. Small energy 
fluctuations for x < 0.8 make it hard to reliably determine Tc max . The asterisks 
represent results from MC Renormalization Group calculations [37| and the solid 
line is the prediction T c (x) — {tanh _1 [e _1 - 45 ^~ Xc ']}~ 1 , converging to the value of 
x c =p c = 0.593 Eg]. 




1 1.5 2 2.5 3 

Temperature [J/k B ] 

Figure 3. The specific heat of the dilute Ising model for different concentrations on 
a L = 22 lattice. For x = 0.75 there is no pronounced peak present 



peaks denoted by £L, corresponding to an "ordered" state (more negative) energy, 
incorporating small clusters, and E + , corresponding to a "disordered" state energy 
incorporating large clustering. A plot of the distribution for the Ising model both 
at T* smg and at T£ w shows clearly sharp single peaks centered approximately at the 

'-'max '-'max ° A tJ A A A ° 

critical energy U c = —\[2J (Fig. [HJ). This supports the uniqueness of the distributions 
in Fig The positions (energies) of the peaks are found to scale with L~ x 39 as seen 
in Fig. and are expected to eventually intersect for a large enough sample. 
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T [J/k fi ] T [J/k B ] 

Figure 4. Calculation of thermodynamic functions for the pure BW model on an 
L = 54 lattice: (a) Internal Energy, (b) Specific Heat, (c) Entropy and (d) Free energy. 
The specific heat displays a very clear pronounced peak at the transition point. 



A comparison between the DOS of the Ising model and the BW model (Fig.IHJ) shows 
a significant difference between the two models. Although they have approximately the 
same number of different energy levels (N — 1 for Ising and N — 3 for BW), the function 
In g, appears to be concave everywhere on the interval [—2, 0] for the Ising model, while, 
for the BW graph this may not be so. This suggests an explanation for the appearance 
of the two peaks which is also consistent with the fact that they have the same height: 
The condition that the distribution will have extrema is satisfied by 

d(\ng)/dE = 1/kT. (4) 

If, at Tc max , Eq. (HJ) has locally, a solution fi(E) = E/kT Cm ^ + G\ tangent to lng at £L 
and E + , and another solution f±(E) = E/kT Cm ^ + C 2 tangent to In g at U C (L) (at the 
shifted critical energy, or the minimum between the peaks), then the distribution will 
have two peaks with equal height satisfying 

p(E_)=p(E + )=e c \ (5) 

as seen in Fig. EJ This is essentially a finite size effect and should be recovered by a 
large enough sample, to give an "Ising like" concave everywhere DOS function, and a 
single peaked distribution as its consequence. 

We further calculated the specific heat for each lattice size and then plotted its 
maximal value C max versus L. We see in Fig. EH a very nice agreement between the 
calculated data and the second order ansatz C max (L) oc L a ^ u , with a/v — 1, even for 
very small lattices (L — 6). 

Another quantity of interest was the so called Binder parameter 001 E] 
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-2 -1.75 -1.5 -1.25 -1 -0.75 



E/N[J] 

Figure 5. Critical distribution calculated at Tc* max for the pure BW model. The 
lattice sizes are denoted by arrows. The L — 120 data suffers from the systematic 
errors resulting from the DOS calculations for large systems. 




-2 -1.75 -1.5 -1.25 -1 

E/N[J] 



Figure 6. Energy distributions at = 2.28948J/fc s , at = 2.27549J/fc B 

and at the exact transition point T c on the same lattice with L = 54. The numbers 
in parenthesis denote: (1) Ising at T§^, (2) BW at T c , (3) fsing at I^f, (4) BW 
at . Note the distribution at the exact transition point (2) with the ratio of 

r ~ 4 between the pronounced peak on the left and the " hump" on the right 26 . The 
asterisk denotes their common critical energy U c = -V2J. 

where (...) stands for the canonical thermal average. When we calculated the Binder 
parameter, (whose plot as a function of temperature is given in Fig.llOj). we saw a sharp 
inverse peak that usually occurs in first order transitions Another manifestation 

of the strong finite size effects is the very precise (though quite unreliable) estimate of 
T c = 2.2696 ± 0.0004 to the transition point we obtained, when performing first order 
finite size scaling theory to the position of B min , Tg m . n . 
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Figure 7. Variation of the energy distribution's two maxima positions with L 1 . The 
"disordered" energies are denoted by (A) and the "ordered" energies by (o). 



Obviously, since the transition is continuous and therefore no ordered and 
disordered states coexist at the transition point, the critical probability distribution 
in the infinite volume limit is expected to be single peaked, causing B min to eventually 
vanish with some exponent and the Binder parameter to take the trivial value of 2/3 
also at the critical point. It was therefore convenient to repeat finite size scaling for 
B min according to 

B mm = ~ - B L- e ^, (7) 

where 9b is an exponent yet to be determined. In Fig. we see the variation of the 
inverse distances t^ 1 ^ = (Tc max — T c ) _1 and tg 1 . = (T BmiM — T c )~ x , correspond to the 
positions of the specific heat maxima and Binder parameter minima, respectively, with 
L. A least square fit gave a slope of 1.529 ± 0.039 for the specific heat temperature and 
1.748 ± 0.025 for the Binder parameter temperature. In accordance with [3T] 

Tew =T C + A L-V» (1 + A\L~ W1 + •••), (8) 

we use the analogy 

T Bmin = T C + B Q L~ x f v (1 + B 1 L- e " + •••), (9) 

where u)\ and 6\ are correction exponents and A ,Ai,B and B 1 are amplitudes 
determined from simulations. It is therefore evident that Ts m . n displays a large 
correction-to-scaling term (0 X ~ 0.25), in contrary to the resulting 1/v scaling from 
the T(7 max fit, which is in fair agreement with the exact 3/2 value, and which is also 
consistent with the scaling of C max . It is also evident, however, from Fig. Upland Eq. ((7|), 
that a and 9b have the same value. Similar exact and simulational calculations of the 
Binder parameter for the Ising model on the same temperature scale (Fig. ITTj) showed 
much broader and less deep minima at , suggesting that these minima vanish with 
an exponent 9b larger than the BW exponent. 
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Figure 8. DOS of BW and Ising models on the L=54 lattice. The function ln g 
appears to be concave "everywhere" in [—2, 0] for the Ising model, while this may not 
be so in the BW case. A plot of In g(E) versus E for a larger BW system (120 x 120) 
is given in the inset. 



250 




Figure 9. Scaling of the specific heat maxima with the length L for the pure BW 
model. The predicted C max (L) oc L behavior is indicated by a solid line. 



5. The dilute BW model 



Let us consider now the ferromagnetic BW model with quenched impurities. The 
Hamiltonian is given by 



n 



-j 



(10) 



We studied systems with lengths L between 18 and 36. We kept concentrations of 
x = 0.8 for L = 18 and of x = 0.9, 0.95 and x = 0.97 for L = 33, fixed, and let them 
vary around x = 0.9 for L = 33. The data for L = 24 was calculated for concentrations 
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Figure 10. The Binder parameter for the pure BW model versus temperature, for 
various lattice sizes, from top (L — 120) to bottom (L — 60) in descending order. 
The Binder parameter is seen in the figure to display an inverse peak whose depth 
decreases as the system size increases. The infinite volume upper bound B^ in was 
estimated using first order scaling theory to B mln (L). 




g i i ~ i i i 

2.2 2.25 2.3 2.35 

Temparature [J/k ] 



Figure 11. Temperature variation of the Binder parameter for the Ising model, from 
top (L = 60) to bottom (L = 24) in descending order. The data in the inset is given for 
the same lattices on a larger scale. The data for L = 54 is calculated using simulated 
DOS. All other data is exact. 

varied around different values from x = 0.85 to x — 0.97. In Fig. El we compare the 
DOS of the pure and dilute BW models. The apparent crossover to a manifestly clear 
second order transition may give rise again to a concave everywhere form of lng, already 
seen for the Ising model in Fig. El The energy levels differ now only in the amount of 
2 J and can take even or odd values for the same lattice size, depending on the vacancy 
distribution. We then performed a calculation similar to that made above for the dilute 



Monte Carlo study of the Pure and Dilute Baxter- Wu model 



12 




Figure 12. Scaling of Te min with the inverse volume of the system for the pure BW 
model. 

Ising model, of T<7 max , to obtain the T c (x) critical line on a lattice with L = 24, and 
then fitted the high concentration data into a continuous (dotted) line (Fig. ITEJ) . All 
the data except for the L = 18 with a vacancy concentration of 0.2, which was, as 
for the dilute Ising model, unreliable because of relatively high dilution, fell very well 
on the dotted line. This may suggest that the critical behavior is rather universal for 
large enough concentrations, because due to the special connectivity of the BW model, 
one would expect smaller energy fluctuations and therefore a larger scatter of data for 
large enough vacancy concentrations, whilst the dilute BW data seems to agree with the 
dilute Ising data for concentrations of x > 0.9. Of course, in order to make definitive 
statements about universality, larger samples would be needed than those used here. 

We performed a rough finite size scaling for the specific heat maxima at a 
concentration of x = 0.9, using the three points measured for L = 33 that were averaged 
and the other data collected for fixed concentrations. Novotny and Landau ^H] predicted 
ajv ~ for a concentration of 0.9. Our results, presented in Fig. IT7I also indicate, at 
least qualitatively, a significant change in a. Since spatial correlations become smaller 
and hence v becomes smaller, the value of a substantially decreases, thus indicating 
an "Ising like" singularity at the finite lattice transition point. Moreover, the Harris 
criterion for the diluted case is hereby confirmed. Another question of interest was the 
influence of vacancies on the nature of the transition. In order to make a statement 
regarding this question we plotted in Fig. EH the energy distribution for different 
concentrations. We see clearly and unsurprisingly that lowering the concentration 
causes the doubly peaked distribution to vanish and become a singled peaked one with a 
narrower width centered away from U c . It may then be plausible to say that in contrast to 
energy fluctuations which become negligible at sufficiently low concentrations, magnetic 
fluctuations increase with increasing dilution and the transition is manifestly second 
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Figure 13. Scaling of the inverse distances Iq and tg., with L, for the BW model. 
The larger slope of the Binder parameter position's fit, may be a result of the large 
correction term d\ (see Eq. . The specific heat data was shifted for ease of reading. 
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Figure 14. Scaling of the quantity (2/3 — -B m i n ) and the specific heat maximum 
Cmax, with L, for the BW model. 



order. 

6. Conclusions 

Our simulations have shown that the WL sampling is a very accurate algorithm. The 
thermodynamic quantities resulting from the calculated g{E), which yield reasonable 
quality critical data, provide good evidence for this. 

Our results show that the pure Baxter- Wu model is strongly influenced by finite 
size effects and corrections to scaling. The scaling of the specific heat maxima is in 
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Figure 15. DOS for a pure (upper curve) and for a dilute (lower curve) BW model 
with x = 0.9, on an L = 36 lattice. 
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Figure 16. Normalized critical temperature for fixed concentrations on different 
lattice sizes for the dilute BW model. 



excellent agreement with the second order form C max oc L a ^ u , even for small lattices, 
and no correction terms are observed. The Binder parameter, however, displays large 
minima for small samples, thus incorrectly could be thought of as a "first order" scaling 
field. It is an "irrelevant" field in the sense that it gives no additional information about 
the universal exponent u, but rather vanishes with an exponent 9b- This exponent is 
also evident in the Ising model and is presumably larger for this model. The vanishing 
inverse peak in both models states that the energy distribution approaches a delta 
function in the thermodynamic limit, although it is essentially non-Gaussian. The 
doubly peaked shape of the latter is rather peculiar. One would usually expect a single 
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Figure 17. Finite size scaling of C max with L for the pure and dilute BW model. The 
data for the dilute model reveals an a exponent close to zero. 




E/N[J] 



Figure 18. Critical energy distribution for various concentrations and critical 
temperatures, calculated on a L — 24 lattice. The numbers in parenthesis denote: 
(1) x = 0.85; T c = 1.74926, (2) x = 0.95; T c = 1.93379, (3) x = 0.97; T c = 2.07671, 
and (4) x = 1 (pure); T c = 2.29164. The distribution is seen to become sharper and 
narrower when the concentration is reduced. 



peaked distribution which becomes narrower, the closer to criticality one is. This shape is 
essentially a finite size effect due to the large fluctuations between the ferromagnetic and 
ferrimagnetic clusters formed in the vicinity of the transition point, and will eventually 
vanish in the thermodynamic limit. The WL method is also very successful when applied 
for the dilute BW model even for small lattices, both in terms of the critical isotherm 
in temperature-concentration plane for a weak dilution, and probability distribution. 
A crossover to a single peaked critical distribution is clearly seen when decreasing the 
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concentration of spins, and a single peaked distribution is evident at a concentration of 
x = 0.85. This is a result of the formation of isolated domains causing relatively small 
energy fluctuations around the critical energy. 

It would be interesting in the future to use larger lattices to confirm our explanations 
of the finite size problems, The relatively high accuracy of the WL method for small 
dilute systems could be applied in the future to study disorder in other models. 
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